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Abstract—Linear embedding via Green’s operators (LEGO) 
is a domain decomposition method in which complex radiation 
and scattering problems are modelled and solved by means of 
interacting electromagnetic “bricks”. We propose an extension of 
LEGO able to handle bodies with anisotropic constitutive param- 
eters and metallic objects (e.g., antennas). Since the anisotropic 
objects are dealt with LEGO, and the metallic parts are treated 
with the electric field integral equation (EFIE), we refer to the 
overall approach as hybrid LEGO-EFIE. The characterization 
of an electromagnetic brick which embeds an anisotropic object 
requires solving a volume integral equation (VIE). Since this 
procedure is carried out for each brick independently, the 
LEGO approach per se is extremely advantageous over the 
direct solution of a global VIE for all the bodies at once. 
Nonetheless, we further mix the hybrid LEGO-EFIE approach 
with the eigencurrents expansion method in order to tackle 
relatively larger problems. The technique is used to analyze 
a reconfigurable plasma antenna array (PAA) comprised of 
magnetized-plasma tubes placed around a two-dipole antenna 
array. 


Index Terms—Antennas, anisotropic media, surface integral 
equations, volume integral equations, equivalence principles, 
Method of Moments, domain decomposition, macro basis func- 
tions, eigencurrents, magnetized plasma, scattering operators, 
plasma antennas. 


I. INTRODUCTION 


Media with anisotropic constitutive parameters exist in 
nature and have found applications such as in the fabrication of 
optical components [1]. Besides, nowadays artificial materials 
can be engineered, for instance, by inserting metallic or 
dielectric inclusions within a host dielectric matrix; then, if 
the distribution and orientation of the underlying inclusions 
is asymmetric, the resulting material is likely to exhibit 
anisotropic behaviour [2], [3]. Also, magnetized plasma dis- 
charges, which have been employed as metamaterials [4] and 
antennas [5], [6], under certain conditions can be described 
macroscopically in the spatial domain by means of a non- 
Hermitian dyadic permittivity [7]. 

Electromagnetic (EM) radiation and scattering problems, 
which simultaneously involve perfectly electrically conducting 
(PEC) objects and penetrable bodies with anisotropic prop- 
erties, are better formulated in terms of coupled volume- 
surface integral equations (VSIE) (e.g., [7]-[11]) because 


Sommerfeld’s radiations conditions [12], [13] are automati- 
cally accounted for by the integral operators. However, it is 
well known that, when solved numerically with the Method 
of Moments (MoM) [14], surface integral equations (SIE) [8], 
[15]-[18] and volume integral equations (VIE) [11], [19]-[24] 
invariably yield dense and possibly large matrices. 

One way of coping with this issue consists of reducing 
the number of degrees of freedom and, thus, the size of the 
algebraic system to be inverted. Said reduction can be realized 
by either breaking the original EM problem into electrically 
“smaller” sub-problems or by adopting specialized basis func- 
tions devised ad hoc for the geometry at hand; combinations 
of both strategies are also possible. An incomplete list of 
the methods that implement these ideas are: the synthetic 
functions expansion (SFX) [25]-[29], the characteristic basis 
function method (CBFM) [30]-[39], and its multilevel version 
[40], [41], the equivalence principle algorithm [42]-[45] and 
the tangential equivalence principle algorithm (T-EPA) [46], 
[47], the generalised surface integral equation (GSIE) method 
[48], the eigencurrents expansion [49], [50], and the linear 
embedding via Green’s operators (LEGO) method [51]-[56]. 

In this paper we propose an extension of the LEGO ap- 
proach to the analysis of EM radiation and scattering problems 
comprised of PEC objects (e.g., antennas) and an aggregate 
of anisotropic bodies, as is sketched in Fig. la. Since the 
anisotropic objects are dealt with LEGO, while the metallic 
parts are treated with the standard electric field integral equa- 
tion (EFIE) so as to gain more generality (the magnetic field 
integral equation alone would, in fact, limit the application to 
closed surfaces) we call the overall approach hybrid LEGO- 
EFIE [57]. 

The basic idea of the LEGO method [54], [56] is the mod- 
elling of a complex and possibly large structure by means of 
simple-shaped EM “bricks”. In this way, the local geometrical 
complexities [55] and even elementary sources confined to 
the inside of the bricks [58] can be treated efficiently, though 
separately, from the interactions which take place among the 
various “distant” parts of a structure. The EM behavior of each 
brick is expressed in terms of surface scattering operators of 
equivalent current densities, whereas the multiple scattering 
that occurs among any two bricks is captured by means of 
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surface transfer operators. By combining these two ingredients 
we can rigorously formulate the original EM problem. For 
this reason, LEGO constitutes a type of domain decomposition 
technique applied to SIEs. 

In principle, there is no limitation to the contents of an 
EM brick, but to handle the aggregate of anisotropic bodies 
of interest here, it is convenient to enclose each body (tagged 
medium ©) into a separate domain, as suggested in Fig. 1b. By 
the same token, the host material (labelled © in Fig. 1b) that 
pads the inside of the bricks can be arbitrary [55]. This feature 
was exploited, for instance, to analyze dielectric waveguides 
in which confinement is provided by electromagnetic band- 
gaps [58], [59]. For the sake of simplicity, though, we restrict 
our discussion to the special case where medium ®© and the 
homogeneous background (labelled ©) possess the same EM 
properties. 

Furthermore, although the EM bricks can have any shape, 
defining bricks with different shapes for a given problem may 
be unwise for computational purposes, or even impractical 
sometimes. Actually, if the bricks have all the same shape — 
but not necessarily the same content — then, the numerical 
calculations of the scattering and transfer operators are facili- 
tated because some intermediate results can be recycled, and 
this helps reduce the computational workload. Secondly, while 
for general distributions of objects in space [60] the shape of 
the bricks is not restricted at all, adopting identical and simple 
shapes may be necessary when one needs to stack the bricks 
tightly in a regular pattern so as to model, e.g., dielectric slabs 
with inclusions [59]. 

The formulation of the antenna problem of Fig. 1 with 
hybrid LEGO-EFIE was outlined in [61] for the special in- 
stance of isotropic bodies, whereas the standard LEGO method 
(i.e., in the absence of PEC antennas) applied to clusters of 
anisotropic objects was described in [62]. In this regard, it will 
be clear further on that the LEGO functional equations take 
on the very same form, irrespective of the nature of the EM 
problem inside the bricks. As a result, the hybrid LEGO-EFIE 
to be discussed here will be based on the same modified EFIE 
we deduced in [61]. By contrast, the calculation of the scatter- 
ing operators — which implies the solution of the scattering 
problem inside the bricks — must perforce be updated in order 
to allow for the dyadic constitutive parameters of medium ©. 
In particular, this could be done by combining the integral 
representation of the EM fields over a brick’s boundary with 
the equivalent variational formulation of Maxwell’s differential 
equations within the brick (e.g., [63], [64]). Alternatively, we 
prefer to adopt an approach based on the MoM solution of 
a VIE for the anisotropic object [62], as this strategy is in 
keeping with the embedding philosophy, which has already 
been developed for PEC and isotropic bodies [54]-[56]. 

It should be noted that, since the anisotropic nature of 
the body inside a brick is rigorously accounted for by the 
relevant scattering operator, which is a surface operator, the 
very application of the LEGO algorithm provides a means for 
reducing the number of unknowns in the context of the MoM 
solution. More precisely, we need to solve at most as many 


independent and relatively “small” VIEs as the number of 
bodies; evidently, this task is far less computationally intensive 
and memory demanding than the inversion of a global VIE 
for the aggregate of objects as a whole. Nevertheless, for EM 
problems modelled with a moderate to large number of bricks, 
the algebraic counterpart of the LEGO functional equations 
can still be too large for inversion with direct solvers, such 
as the LU factorization [65]. Therefore, in an attempt to 
overcome the memory limitations that large matrices would 
pose, we combine the hybrid LEGO-EFIE approach with the 
eigencurrents expansion method (EEM) [54], [57], [61], [66]. 

Simply put, in the EEM the eigenfunctions of a brick’s 
scattering operator are employed as a set of locally entire- 
domain basis functions to expand the unknowns (i.e., equiva- 
lent surface current densities) introduced over the boundary of 
a brick. Since ordinarily only few lower-order eigenfunctions 
[57] are sufficient to accurately represent the unknowns, the 
EEM allows compressing the algebraic system of LEGO and 
hence, speeding up the matrix filling phase and the inversion. 

The rest of the paper is organized as follows. In Section 
II, we review the basic concepts of LEGO (II-A), introduce 
the relevant unknown current densities through the EM equiv- 
alence principles (II-B), derive the LEGO equations for the 
bricks (II-C), formulate the EFIE for the conducting parts of 
the problem (II-D), and deduce the modified EFIE (I-E). In 
Section III the derivation of scattering and transfer operators is 
discussed with a focus on bodies with anisotropic constitutive 
parameters. The numerical solution with the MoM and sub- 
sectional basis functions is outlined in Section IV, whereas the 
compression of the algebraic counterpart of the LEGO equa- 
tions with EEM is the subject of Section V. The validation of 
the numerical implementation and the convergence of the EEM 
are presented in Section VI. As an example of application of 
hybrid LEGO-EFIE to a complex EM problem, in Section 
VII we investigate the radiation properties of a reconfigurable 
plasma antenna comprised of a two-dipole linear antenna array 
surrounded by cylindrical magnetized-plasma tubes. 

Finally, a time-dependence in the form exp(jwt) for fields 
and sources is assumed and suppressed throughout. 


II. FORMULATION 
A. Description of the problem and LEGO model 


We are concerned with the solution of an EM radiation 
problem (Fig. la) which consists of a metallic multi-port 
antenna and Np > 1 identical penetrable bodies immersed in 
a homogeneous unbounded background (medium ©). Our goal 
is to determine the network parameters (e.g., the impedance 
matrix) and the radiation pattern of the antenna-objects system. 

The objects — which can have arbitrary shape and position 
in space — are made of an anisotropic material (medium ®) 
endowed with dyadic permittivity (€3) and permeability (723); 
let Vo denote the region of space occupied by an object. As is 
customary, we assume the antenna and the other metallic parts, 
if any are present, to be perfect electrical conductors (PEC). 
We indicate with Va the volume occupied by the conductors 
and with Sa = OV, the boundary thereof. Besides, for the 
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Fig. 1: EM problem and LEGO bricks: (a) a metallic antenna operating in 
the presence of Np anisotropic bodies (medium ©); (b) in accordance with 
the LEGO approach, each object is embedded inside an EM brick Dz, k = 
1,..., Np, in turn described by means of a scattering operator Sgk. In this 
work we examine the special instance of background material (©) and host 
medium (@) endowed with the same constitutive parameters. 


sake of simplicity we model the antenna excitation by means 
of the delta (or voltage) gap generator [67], [68]. As a result, 
since we set the voltage at each antenna port and we compute 
the current flowing therein, the natural network quantity, which 
can be obtained by solving the EM problem, is the admittance 
matrix. 

We are now ready to apply LEGO to the EM problem of 
Fig. la. We start by enclosing the Np anisotropic bodies in 
as many simple-shaped bricks Dk, k = 1,...,Np, and we 
denote the boundary of Dy with Dp. The unit normal fi, to 
OD, is taken as pointing inwards Dy (Fig. 1b); correspond- 
ingly, the positive (negative) side of OD, is indicated by OD{ 
(OD,, ). However, here the distinction between ôD} and 0D, 
is immaterial, since OD% constitutes a mathematical separation 
surface [54] rather than a material interface (cf. [55]), and the 
EM properties of medium © and © are the same by hypothesis. 

We suppose that D; has been characterized by means of 
the scattering operator Skk, whereas the multiple scattering 
among Dz, Dy, n = 1,...,Np, n Æ k has been accounted 
for through the pair of transfer operators Tkn and Tng [54]. To 
keep the exposition lucid, though, we postpone the derivation 
of Skk, Tkn, and Tng until Section TI. 

As we are concerned with a cluster of Np identical bodies, 
all the scattering operators can be made equal to each other, 
provided we also consider bricks with the same shape, as we 
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Fig. 2: Equivalent EM problems: (a) The SEP is repeatedly applied on a 
surface that tightly wraps one brick at a time (say, Dz) so as to preserve 
equivalence inside that brick and to remove the antenna (Va) as well as the 
remaining bricks (Unx DPn); (b) The SEP is applied on a disjoint surface that 
tightly encloses all the bricks and the antenna so as to preserve equivalence 
outside the bricks and the antenna region, and to remove the antenna (Va) 
and the bricks (Uz Dx). 


have argued in the Introduction. Nevertheless, the equations to 
be derived in the following are completely general, inasmuch 
as they depend neither on the shape nor on the content of Dx. 


B. Equivalent problems and unknown current densities 


The theoretical framework for the definition of the EM 
bricks and the formulation of the EM problem of Fig. la is 
the Surface Equivalence Principle (SEP) [14], [69]. Thanks to 
the SEP, equivalent problems in a given region of space can 
be built as long as suitable surface electric (J) and magnetic 
(M) current densities are placed on the separation surface. 

For instance, we invoke the SEP Np times on OD; k = 
1,..., Np, as illustrated in Fig. 2a, in order to obtain Np new 
problems that are equivalent to the original one inside D,, 
whereas the outside of Dy is the excluded region. With each 
step we end up introducing a set of equivalent total incident 
surface current densities, namely, 


i — a i 
k,tot = Nk X Hi tot» 


(la) 


k,tot = Ek,tot X De; (1b) 


where Ei. sot Hi yop represent the total (unknown) incident 
EM field on the boundary OD x, and fi; is the inward pointing 
unit normal thereon. Moreover, the above currents radiate in 
the presence of the object inside Dk only, provided we “fill” 
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the excluded region with a material that has the same EM 
properties as those of medium ©. 

A further application of the SEP, simultaneously on a 
disjoint surface comprised of the boundaries we 2, OD, plus 
the antenna surface Sa (see Fig. 2b), allows us to define an 
equivalent problem in the region of space outside Oe Dr 
and V4. In the process we introduce i) Np sets of equivalent 
scattered surface current densities, viz., 


(2a) 
(2b) 


where E$, H; is the (unknown) scattered EM field on 0D,, 
and ii) the equivalent surface current densities 


Ja = ña x Hy, 
Ma = Ea x fig = 0, 


(3a) 
(3b) 


where Ea, Hy are the total electric and magnetic fields on 
the antenna surface Sa, and ña denotes the unit normal 
thereon. Lastly, upon replacing the inside of the bricks and the 
conductors with a material endowed with parameters £1, 41, 
we are left with J, M;, and Ja radiating in a homogeneous 
unbounded medium. 

We observe that J ` ae Mi, tot (Fig. 2a) radiate the incident 
EM field due to all external sources towards the inside of 
D,. For the EM problem of Fig. 1b, such sources will be 
both (2a) and (2b) flowing on the boundaries of the remaining 
Np-—1 bricks and (3a) induced on the antenna conductors Sa 
(Fig. 2b). Conversely, each pair J;, Mè produces the scattered 
EM fields that are actually radiated by secondary (polarization 
or magnetization) volume currents flowing in the region Vo 
occupied by the anisotropic object within Dķ. 


C. LEGO equations for the bricks 


The relationship between Jise Mi, tot and J}, M}, just 
defined, is provided by the scattering operator Skk, as follows: 
dk = Skk aktov On Dk, (4) 


where gi, sot Gj, are the abstract 2 x 1 vectors 


i Vndi, tot | 
dk,tot = ae ; (Sa) 
pace | =M} tot/v nı 
s [| vmig | 
qk = g 7 (5b) 
i | —M;//m 


and 7 = y uı/E£1 is the intrinsic impedance of medium 
©. The latter normalization factor endows Oh knee qy With the 
physical dimension of a power wave, i.e., w!⁄2/m. In light of 
(4), Skk can be rightfully regarded as a generalization of the 
scattering matrix in network or waveguide theory [70], [71]. 

To proceed, we split the total equivalent incident currents 
into the contributions of the various sources, namely, 


. . No : 
dktot = Te + > Fein): (6) 
er 


where Gs represents the equivalent incident currents due to J A 
on the antenna conductors, and the nth term of the summation 
constitutes the equivalent incident currents due to q), flowing 
on the boundary of D,,. The link between equivalent scattered 
currents on OD,, and additional incident currents gi, (n) On OD; 
reads 


on Dk, nÆ k, (7) 


where T;,, is the transfer operator between OD,, and OD, (see 
Section II). 

Next, upon substituting (7) into (6) and the latter into (4), 
and then formally solving for qi, and rearranging, we arrive at 
the set of Np coupled functional equations 


k(n) ae Tkn In ; 


Np , 
Siete — >, Tend, = Ges 


n=1 
NÆR 


k=1,..., Np, @®) 


which can be cast in a compact abstract matrix form as: 
Stg = é, on ôD: U... U Dyp, (9) 


where q% are Np x 1 abstract vectors 


qi 
TAAIE NE (10) 
INp 
and 
S~! = diag{S;} } — T, (11) 


is the total inverse scattering operator of the aggregate of 
the bricks [54]. Additionally, the total transfer operator T, an 
abstract Np x Np matrix, is defined as 


0 -T2 itis 
=T 0 —T2Np 
T= : , i ; (12) 
—Tnpi —Tnp2 °°: 0 


and, since Tkn Æ Tnx, it is not symmetric (Section IM). 

The total inverse scattering operator in (11) rigorously and 
elegantly captures the multiple scattering that occurs among 
the Np anisotropic bodies. And yet, (11) is just a formal 
result, because first of all Sz, and its inverse cannot be 
obtained analytically for general geometries; secondly, Sa. 
may not be defined at all when the scattering operator happens 
to be singular. Nevertheless, we shall show in Section V 
that the numerical solution of (9) is by no means critical, 
provided we apply the eigencurrents expansion and organize 
the calculations properly. 


D. Electric field integral equation for the antenna 
We now turn our attention to the antenna part of the EM 
problem. By enforcing the boundary condition (3b) for the 
total tangential electric field on Sa we arrive at the electric 
field integral equation (EFIE) 
Np 
E$ +E) +) Eå 
k=1 


=0, on Sa, (13) 








tan 
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where EÑ is the impressed field supplemented by the external 
generator in the delta-gap approximation [67], [68], E% is the 
secondary field produced by Ja on Sa, and Ej, represents 
the field radiated by the scattered currents gj on OD,;. The 
subscript “tan” implies that (13) must hold true for the field 
components that are tangential to Sa. 

Since the voltage gap generator [67] consists of the quasi- 
static approximation of the electric field existing at the open 
end of the coaxial cables connected to the antenna, E$ is 
effectively non-zero only at the antenna ports. In fact, the 
voltage gap model can be applied successfully as long as the 
actual antenna port is small as compared to a) the operational 
wavelength in the background medium and b) the antenna size. 
As an example, in the case of a single-port antenna 


[EX lean = VaO. on Sa, (14) 


tan YA? 


where VG is the intensity of the voltage generator, y4 is a line 
on Sa that defines the antenna port, Ô is a unit vector tangen- 
tial to Sa and locally perpendicular to y4, and 6,, denotes 
a one-dimensional Dirac distribution on y4. Furthermore, in 
light of (14), there is no direct contribution from the generator 
onto the aggregate of bricks; this is reflected in the absence 
of a corresponding equivalent incident current in (9). 

The tangential part of the electric field produced by Ja on 
itself is given by 





[Eà lan = mLanrJa, on Sa, (15) 
where 
Lawr{o} = ih f dr'a (R)I {o} 
Sa 
JVs 2 , 
= d*r'Gy(R)V,-{o}, resa, (16) 
ky Sa 


is the standard EFIE operator [14]. Besides, kı = w,/ey p41 
represents the wavenumber! in medium ©, R = |r — r'| is the 
distance between field (r) and source (r’) points, and 


exp (—jkiR) 
= —— 17 
Gy (R) AnR , (17) 
is the pertinent 3-D scalar Green function. Finally, I, = I — 


Haha is the surface unit dyadic tangential to Sa, Vs = I.-V, 
and V| = —V, means differentiation with respect to r’. 


E. Formulation with a modified EFIE 


To finalize the formulation of the EM problem we need 
to include the multiple scattering phenomenon which occurs 
between the antenna and each one of the Np bodies. 

For instance, we observe that the antenna current J4 radi- 
ates the EM field Ej, Hi on the boundary of Dk. This can 
be written symbolically as 


. 0 
F} = ; = PraJd 
tk | Vm Hi, | Vm kAYA, 


where the surface integro-differential operator Pa is called 
the propagator from Sa to OD;, (see Appendix A-C); besides, 


(18) 


!The wavenumber kı should not be confused with the brick index k. 


the subscript ‘t’ stands for “tangential” to OD;. To convert the 
fields impinging on OD, into the equivalent incident currents 
in (6), we resort to another propagator, as follows: 


kk Te = Fik, on Dx, (19) 
where Phi is a 2 x 2 abstract matrix of integro-differential 
operators on OD;, which is given explicitly in Appendix A-A. 

Note that a null vector appears in the definition of Fi, 
above, because the propagator Ppi; has been conveniently 
defined in such a way that only the incident magnetic field 
is required to produce the corresponding incident currents. 
Thanks to this expedient step followed, in the context of the 
numerical solution through the MoM and a symmetric inner 
product, the algebraic counterpart of Piy happens to be a 
symmetric matrix. It is also worthwhile mentioning that the 
calculation of É, through (19) is preferable over the definition 
in terms of twisted tangential fields, namely, 


J; 


Mi = Ei x fig, 


fi, x Hi, (20a) 


(20b) 


since the latter equations are likely to yield less accurate results 
when evaluated numerically. 
Now, solving (18) and (19) for di yields 


d, = vm (Pip) Prada, 


which mathematically accounts for the interaction between the 
antenna and the brick Dy. 

The equivalent scattered currents q}, flowing on OD, as a 
result of the presence of the body inside Dx, in turn produce 
the EM fields Eñ, H4,, on the antenna surface. Since only 
the electric field enters (13), we can express Eį, symbolically 
as 


on Dk, (21) 


[EA klhtan = V nı Pak Ge 


where the surface operator Pa, is called the propagator from 
Dp to the antenna, and it is a 1 x 2 abstract matrix of integro- 
differential operators on OD,, as detailed in Appendix A-C. 

With the aid of (10) and (22) we can cast the combined 
contribution of the Np bricks onto the antenna in a compact 
form as follows: 


on Sa, (22) 


Np 
Yo TEXaltn = Vii Paog’, on Sa, (23) 
k=1 
where Pao is a 1 x Np abstract matrix with entries 
(Pao)k = Pag, (24) 


and qê is given by (10). By inserting (21) into (9) and solving 
for që we obtain 


€ = vm SToaJa, 


where Toa is an Np x 1 abstract matrix of operators with 
entries 


on 0D, U... U ƏDyp, (25) 


(Toa)& = (Phx) Pra. (26) 
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Finally, inserting qê as given by (25) into (23), and the result 
into (13), yields the modified EFIE for the current density Ja, 
viZ., 


m (Lant + Leco) Ja =—[E§]tan, on Sa, (27) 


where the operator 


Linco = PaoSToa, (28) 





factors in the effect of the aggregate of bricks (and hence, of 
anisotropic bodies) on the antenna current J4. In this regard, 
if we were able to cast Lant + Luego in the form of a 
radiation integral over S4 and to single out the corresponding 
kernel explicitly, we would have obtained the Green’s function 
relevant to the electric field produced by an elementary electric 
current density in the presence of the Np objects. 

The fact that (27) and (28) have exactly the same functional 
form as [61, Eq. (5)] — which were derived for antennas and a 
cluster of isotropic objects — is a welcome consequence of the 
modularity of LEGO. Indeed, it is evident that the assumption 
of anisotropic objects has an impact solely on the calculation 
of the scattering operators Skp. 

Once the antenna current density Ja has been computed 
from (27), we can determine the equivalent scattered densities 
q from (25). Both Ja and q° contribute to the radiated EM 
field. 





III. SCATTERING AND TRANSFER OPERATORS 


In this section we derive explicit, though formal, expressions 
for Skk and Tkn, Tnk in terms of suitable integral operators. 
We recall that in deriving S;,, we assume that media © and 
®© possess the same EM properties. To gain more generality, 
we suppose that the object inside D; has dyadic permittivity 
é€3 and permeability m3. Also, by invoking linearity and 
superposition we can work with a solitary brick D; in medium 
© illuminated by incident EM fields Ei, Hi, which are 
generated by some external source of radiation. In fact, we 
may as well assume that the equivalent incident currents q 
have been obtained by inverting (19). 

To begin with, we observe that the EM field produced by 
g, in the region Vo occupied by the object inside D; can be 
written formally as (Fig. 3a) 


F= | D,/(e1 vm) 
° L Boy /H1 


with Dİ, B', the incident electric and magnetic flux densities 
within Vo. The propagator Pox is a 2 x 2 matrix of dyadic 
surface operators whose kernel is the Green’s function of 
media © and © (see Appendix A-B). 

Next, we apply the Volume Equivalence Principle (VEP) 
and we replace the object with the equivalent electric and 
magnetic volume current densities (e.g., [22]) 


Jo = jw[I— (Es()/e1)""] Do 


| = Pordi, in Vo, (29) 


=jw@elr)- Do, rE Vo, (30a) 
Ms = jw[I- (f3(r)/11)~*]- Bo 
= jw@p(r)- Bo, rE Vo, (30b) 


where Œe and @» are the contrast dyadics, and D, and Bo 
are the total electric and magnetic flux densities, respectively, 
within Vo. Accordingly, Do, Bo can be obtained from Fi by 
solving the VIE (Fig. 3b) 


Xoo qo = Fi, in Vo, (31) 
where 
Do 
=| eh (32) 
ym Bo/ Ha 


and Xo, is a 2 x 2 matrix of volume operators over Vo (see 
Appendix B). The volume currents Jo, Mo in turn generate 
a scattered EM field on the Dp. The latter can be expressed 
mathematically in terms of Do, Bo as (Fig. 3c) 


Ss a 


tk = | = Pkoqo, on Dk, (33) 


| 0 
vm Hir 
where the propagator Pko is a 2 x 2 abstract matrix of dyadic 
volume operators involving the Green’s function of media © 
and ®© (see Appendix A-B). The superscript ‘s’ (short for 
“scattered”) signifies that the fields in question propagate away 
from the object and out of the brick. 

Eventually, a further application of the SEP over OD; 
enables us to replace the scattered field Fẹ with equivalent 
scattered densities (Fig. 3d), as follows: 


Pind, = Fh, on OD, (34) 


where qj, is the abstract 2 x 1 vector introduced in (5b). The 
propagator P’, — which plays a similar role to Pİ, in (19) 
— is defined so as to require only the magnetic field HẸ, to 
determine q}, and this accounts for the presence of the null 
vector in F$, in (33). In this way, as already observed for P;,,, 
the algebraic counterpart of P;,, is also a symmetric matrix. 

Now, by formally solving (29), (31), (33) and (34) for qj, 
as a function of qd we find that? 


di. = (Pix) 1Pko(Xoo) Par Gk = Skx Ges 


which provides us with a suitable recipe for computing the 
scattering operator Skk. The propagators Pko, Pok and the 
integral operator Xo. are given in Appendices A-B, B for the 
special instances of either a dielectric body (with m3 = pI) 
or a magnetic one (with €3 = cD). 

The interaction between bricks Dk, D,, can be obtained as 
follows. In a similar fashion to (19), we turn the incident fields 
onto OD; (or OD,,) into equivalent surface current densities 
dein) (or dni k)) Symbolically, we write 


(35) 


n £k, 


kÆ£n, 
where the propagators Pi, Pi,,, are the same as in (19). 
Besides, the propagators Pz,, Pnę (given in Appendix A-A) 


are 2 x 2 abstract matrices whose non-null entries are dyadic 


(36a) 
(36b) 


ae de(n) = Pkn Dns on OD,, 
Pin In(k) = Pak Teo on OD,,, 


This expression for Sgp differs from the one in [54, Eq. (11)] for a minus 
sign, which is due to different definitions adopted for the operator Xoo. 
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Fig. 3: Derivation of the scattering operator of Dy in the case when host (®©) and background (©) medium have the same EM properties: (a) The equivalent 
incident currents qj, on ODx produce incident fields FG in the region occupied by the object; (b) The object “reacts” with polarization and/or magnetization 


currents qo; (c) The currents go produce scattered fields F§, 
integro-differential operators. The formal inversion of (36a), 
(36b) yields 


dalk) = (Pin) Pan dk = Tnk dk» 


n#k, 
k#n, 


where the dimensionless transfer operators Tkn, Tnk map 
scattered currents on OD,, ;, onto incident currents on OD; ». 
The propagators involved in (37a) and (37b) depend solely 
on the shape of a brick. Therefore, regardless of the content of 
a brick, if D;, and D,, have the same shape, then Pi, = Pin, 
hence, we need only compute either one. By contrast, even 
in the case when OD; = Dp, it turns out that Pen Æ Pre. 
For these reasons — but also because Pi, and Pi, do not 
commute with P;,, and Png — no simple relationship exists 
between Tkn and Tng, and the total transfer operator T in (12) 
is not symmetric, as we anticipated. It can be shown that the 
algebraic counterparts of the non-null entries of P;,, and Prk 
(see Appendix A-A) are transposes of each other, provided the 
MoM is implemented by using a symmetric inner product. 


(37a) 
(37b) 


IV. NUMERICAL SOLUTION WITH MOMENTS METHOD 


We now describe the application of the baseline Method of 
Moments (MoM) in the Galerkin’s form [14] for the calcu- 
lation of scattering and transfer operators in (35) and (37a), 
(37b), the object-to-antenna propagators (24) and the antenna- 
to-object transfer operators (26). Furthermore, we derive a 
weak form corresponding to (25) and (27). 

We begin by modelling the surface of the bricks and the 
conducting elements with 3-D triangular tessellations, and the 
anisotropic objects with tetrahedral meshes. Then, we express 
the surface electric and magnetic current densities flowing on 
the boundary of a brick as [54] 


i — fkp(r) y/n T Jip 
dk = 2 | grplr)M i AVT | , r €8Dk, (38a) 
sS io fkpr) y/n iJ; p 
q; = ` i P za (ME, [am | , reðôD,, (88b) 


on Dx; (d) The scattered fields are replaced with equivalent scattered currents h- 


where {f O hE and {gr pE constitute two sets of Np 
sub-sectional real divergence-conforming surface vector basis 
functions [72] associated with all the edges of the underlying 
3-D triangular-facet mesh. Here, we choose frp = Skp, al- 
though the functions used to express the electric and magnetic 
current densities need not be the same. 

Likewise, we expand the flux densities within an anisotropic 


body as [62] 


S| 


m=1 


V(t) Dm/(E1/™) 
Vm(r)/ Mm Bm/H1 Í 


rE Vo, (39) 


where {vm (r)}¥2] represents a set of No sub-sectional real 
divergence-conforming volume vector basis functions [73]. 
In particular, if the body is either electric or magnetic, (39) 
will reduce to an expansion for just Do or Bo, respectively 
(see Appendices A-B and B). To account for the fact that 
the normal components of D, and B, do not vanish on 
the boundary OVo of the object, the functions vm are to be 
associated with all the triangular facets of the tetrahedral mesh. 
Finally, we express Ja on Sa as [61] 


Na 
Jar) = So hi(r)h, re Sa, (40) 
l=1 


where {h;(r)} 4 are N4 real divergence-conforming surface 
vector basis functions [72] associated with the inner edges of 
the relevant triangular tessellation. 

To finalize the application of the MoM we also need three 
L? symmetric inner products defined as 


(a, F)ap, = dr a: F, ac {fkp; Birhan (41a) 
OD; 
(a,F)y = | @ra-F, ac {vm}%2,, (41b) 
Vo 
(a,F)s,= | d’ra-F, ac {h)}%4, (41c) 


Sa 


with F a complex vector field in Dk, Vo and Sa, in (41a), 
(41b) and (41c), respectively. 
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For ease of manipulation, it is convenient to collect the 
expansion coefficients in (38a), (38b), (39) and (40) into 
column vectors [q,], [g$], [do], [Ja], as follows 


Fv Di/(e1/m) 


Jye VM Dno/(&1yM) 


RS l + a= , (42) 
[a | _ Si) AR [qo] Jm Ba /m 
_ gi, Vm VmBno/ HA 
qi 
[Ja] = (43) 
In, 


Following the manipulations above, we proceed to deter- 
mine the weak form of equations and operators. 

For instance, we insert (38a) into (29), (39) into (31), and 
we project them onto Vo with (41b); then, we insert (38b) into 
(34), and we project it onto OD; with (41a). Upon solving the 
trio of resulting equations for [qj] in terms of [gi], we arrive 
at the expression for the algebraic scattering operator 


[Sie] = Peel [Pro] [X00] * [Por], (44) 


where with transparent notation each matrix represents the 
algebraic counterpart of the corresponding integral operator 
in (35). As regards the size of the matrices in (44), [P] 
is 2Np x 2Np, [Xoo] is 2No x 2No, and [Pro], [Pox] 
are 2Np x 2No, 2No x 2Np, respectively. However, if the 
object within Dx is either electric or magnetic, then [Xoo] is 
No x No, and [Pro], [Pox] are 2Ne x No, No x2Np, because 
either M, or J, is zero in Vo. 

By inserting (38b) and a similar expansion for dh, a into 
(36a), projecting onto OD; with (41a), and formally solving 
for [q}] versus [dkny we can express the algebraic transfer 
operator as 

[Ten] = [Peel [Pin] > (45) 
where both matrices are 2N p x 2N p; the expression for [ Tnk] 
ensues from (45) by swapping the indices n, k. 

By substituting (38a) into (19) and (40) into (18), projecting 
onto ôD, with (41a), and solving for [qi] with respect to [Ja] 
we get 


[Toar] = [Pi] t [Pra], (46) 


where [Ppa] is size 2Np x N4. Similarly, inserting (38b) into 
(22) and projecting onto Sa with (41c) yields [Pap], whose 
size is N4 x 2Np. The matrices [Toax]| and [Pax] are needed 
to form [LyEcol, the algebraic counterpart of (28). 

By inserting (40) into (15), and projecting onto Sa with 
(41c), we obtain [Lan], the algebraic counterpart of the EFIE 
operator (16). More specifically, the source term corresponding 
with (14) and the right-hand side of (27) is a column vector 


[EÑ] with N4 entries given by 
ER, = | rh(r)- E$ (r) 
Sa 


= Ve | dr hj(r)- Ô, (47) 
YA 

which, as a result of the delta-gap generator model (14), are 
non-null only for those test functions among h;(r) that are 
associated with the inner edges spanning the line ya in the 
triangular-facet mesh of Sa. 

Now, by making use of (44) and (45) the algebraic total 
inverse scattering operator (of size 2Np Np x 2Np Np) can 
be written as 


[Su] [T12] — [Tinb] 
isj-? = = [Fa] a = [Taw (48) 
a [Tvp1] = [Tn p2] [SNoNp] ` 


whereas the algebraic equivalent of the LEGO operator (28) 
reads 


[LrLeco] = [Pao] [S] [Toa], (49) 
with 
[Pao] = [[Pai]--- [Pan]; (50) 
[Toa1] 
[Toa] = : : (51) 
[Toan] 


and [S] is the inverse of [S]~’. 
Lastly, in light of (47), (49), (50) and (51), the weak form 
of (25), (27) reads 





[°] = vm [S] [Toa] [Ja], (52) 
m ([Lant] + [Lreco]) [Ja] = — [E3], (53) 


where [q] is a column vector with Np block entries [gj,] as 
defined in (42). 

As long as the number Ny of basis functions h;(r) is not 
too large, we can solve the system (53) (of rank N4) by 
means of direct solvers, such as the LU factorization [65]. The 
real challenge, though, is posed by the very determination of 
the LEGO matrix [Leco], because that computation entails 
an inversion of [S] in (48), as can be gleaned from (49). 
In fact, even if we suppose the scattering matrices [Skp] are 
invertible so that [S]~" can actually be built explicitly, [S]~* 
may already become a large matrix, even for a moderate 
number Np of EM bricks in the model. Worse still, if [Spp] 
is singular, then the representation (48) — being just formal 
— cannot be used to determine [S] in the first place. These 
considerations also apply to (52) which we need to retrieve 
the coefficients of the scattered currents on UN? ODE Thus, 
to cope with the occurrence of a large and possibly undefined 
matrix ISI}, we adopt the eigencurrents expansion method 
(EEM). 
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V. EIGENCURRENTS EXPANSION AND COMPRESSION 


The application of the eigencurrents expansion method 
(EEM) to the solution of the LEGO equation (9) alone was 
extensively described in [54], [66]. Here, we follow the same 
strategy to compute the matrix product [S] [Toa], which enters 
the algebraic LEGO operator (49) and the expression of the 
scattered currents coefficients [q°] in (52). Therefore, the main 
goal is the efficient and stable calculation of [S]. 

To begin, we consider the matrix of 2N eigenvectors [Vix] 
of the algebraic scattering operator (44) defined through 


[Six] [Vix] = [Vex] diag {A}, 


where aP, q= 1,...,2Np, denote the eigenvalues of [Skk]. 
We also assume that the eigenvalues are sorted in descending 
order, i.e., [vs > Beal and that this ordering is reflected 
on the columns of [V;,] in an obvious way. Upon combining 
the 2p elements of each eigenvector with the basis functions 
fp, Skp» as is prescribed by (38a), we obtain on OD, a set of 
surface current densities which we refer to as the eigencurrents 
of Dp, in light of (54) and the role played by [Skk]. 

In the eigencurrents expansion method (EEM) we employ 
the eigencurrents as local entire-domain basis functions on 
Dp. The desired compression of [S]~* is then achieved 
because only the eigencurrents associated with the larger 
eigenvalues aP (usually only a few) must be retained to 
expand the currents on 0ODx,. 

This is accomplished by first expressing [S] 7} in the basis 
of the eigenvectors of [Skk], as follows: 


[57 = [Vv] [S1 [v], 


(54) 


(55) 


where [V] is a block-diagonal matrix with Np diagonal blocks 
given by [Vk]. More specifically, in view of (48) and (54), 
we can write 


(l en = eo > o 


— [Ver] " [Ten] [Van], k £n. 


Now, [S]~! would be an exactly diagonal matrix, if we had 
used its true eigenvectors, but the latter are just the quantities 
we cannot compute either because isi is too large, or 
because [Skk] is singular, or both. As it turns out, though, 
most of the entries of (55) are relatively small and can be 
safely set to zero. 

To this purpose, we observe that the eigencurrents of Dk 
can be separated into two complementary sub-sets: coupled 
and uncoupled. The former — corresponding with the larger 
eigenvalues aP, q= 1,..., Nc — produce substantial radi- 
ation, and hence they contribute to the multiple scattering that 
takes place among the bricks. On the contrary, the uncoupled 
eigencurrents — associated with the smaller and possibly null 
eigenvalues ag, q= Nc +1,...,2Np — generate scattered 
fields that remain mostly confined around Dx; hence, they are 
only necessary to describe the scattered currents on the very 
same brick. 

To proceed with the inversion of [5] 
and columns as follows: 


= we rearrange rows 






= =i = ami 
[Scc] [Avu] 
NEGLIGIBLE! RECIPROCAL 
CONTRIBUTION $< ENTRIES || . HIGHER-ORDER 
OF THE COUPLED \ EIGENVALUES 
EIGENCURRENTS ASSOCIATED WITH 
THE UNCOUPLED 
EIGENCURRENTS 


Weel! (Sresceeseloer esse eseeueeees. 


Fig. 4: To illustrate the partition and subsequent compression of the total 
inverse scattering operator [s]7t in the basis of the eigencurrents. 


1) The entries relevant to two coupled eigencurrents are 
shifted to the upper left corner of the matrix; 

2) The entries corresponding with two uncoupled eigen- 
currents are displaced to the lower right corner of the 
matrix; 

The entries pertinent to a pair of coupled and uncoupled 
eigencurrents are moved either to the upper right part or 
to the lower left part of the matrix. 


3 


wm 


Symbolically, we can write 


[Soc] [Scu]! 


[Suc]! [Suu]! 


where the subscript ‘C’ (‘U’) stands for coupled (uncoupled) 
and [|P] represents a suitable 2NeNp x 2NrNp permuta- 
tion matrix [65]. In accordance with the notion of coupled- 


uncoupled eigencurrents we observe that: 
7 


[Â = [P]" [5] [P] = 





| » (57) 


1) the matrix [Soc is dominant; 

2) the off-diagonal matrices [Scy]~! and [Syc]~! must be 
relatively small; 

the off-diagonal entries of [Syu]7 
small or null. 


3 1 


wm 


must be relatively 


Based on these observations, we approximate (57) as 


«1 fie wo 
A =| 0 ae si 


with [Auu] representing a diagonal matrix that contains all the 
eigenvalues associated with the uncoupled eigencurrents; the 
reduction process is graphically illustrated in Fig. 4. In this 
way, the calculation of [S] has been reduced to the separate 
inversion of [Scc]~! and [Ayu]7!. 

The matrix [Soc]! has the rank Nc Np, which is relatively 
small as compared to the rank of the original uncompressed 
inverse scattering operator (48); thus, we can obtain [Scc] with 
the LU factorization [65]. Besides, the diagonal of [Scc]! 
is comprised of the reciprocal eigenvalues germane to the 
coupled eigencurrents, as is evident from (56). Since these 
eigenvalues are relatively large — although rvs” <1— 
their numerical inversion is stable, and hence, [Soc]! is 
well defined. In addition, we need not build [Ayu]~! at 
all, since from (55) we already know the eigenvalues of 
[Skk]. Therefore, the occurrence of very small or possibly null 
eigenvalues does not pose any numerical issue. 
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To determine the algebraic LEGO operator (49) we first ex- 
press the matrices [Pao] and [Toa] in the basis of eigenvectors 
of [Skk] as well and then shuffle the entries thereof with the 
permutation matrix [P]. The result reads 


[LiEco] = 
= [Pao] [V] [P] [$] [P]" [V+ [Toa] 


=[[Pao.cl [Pao.ul] [S] ee 


~ ([Pao,cl[Scc][Toa,c] + [Pao,u][Avu][foaul, 69 


where we have made use of (58) in computing [S]. The second 
contribution in the rightmost r.h.s. of (59) can be neglected 
altogether, because the uncoupled eigenvalues are relatively 
small. From a physical standpoint, the latter implies that the 
uncoupled eigencurrents generate fields that do not propagate 
towards the antenna. 

For the calculation of [q°] in (52) we proceed in a similar 
fashion, namely, 


l°] = vm [V] [P] [$] [P] [V] [Toa] [Jal 


= vin VIPIS |AS ua 
[Toa,u] 
a [Sco] [Toa,c] 
= vm [V] [P] Aca] Toast [Ja] , (60) 


where again the matrix (Auu][Toa.u] can be neglected. 

If all the scattering operators are equal to each other, then 
we have to perform the spectral decomposition (54) and the 
formal inversion of [Vp] only once. Besides, the total number 
of coupled eigencurrents will be Nc Np, which is the size of 
[Seo]74, the only matrix we really have to invert to compute 
[LiEco] and, at a later stage, the scattered current coefficients 
[a°]. 

We conclude this section by observing that, from a compu- 
tational viewpoint, it is desirable to use as few eigencurrents 
as possible, so as to keep the size NcNp of [Soc]! at 
bay. Obviously, the number of coupled eigencurrents cannot 
in any case exceed the total number of basis functions 2N p 
introduced on OD; to expand the equivalent current densities 
thereon. Furthermore, the very form of (44) suggests that [Sk] 
may be rank deficient, if the size of [Xoo] is smaller than the 
size of [P$]. Therefore, the maximum number of coupled 
eigencurrents is 


Nc,max = min{2Npr,2No}, (61) 
if the embedded object is both magnetic and dielectric or else, 
Ne,max = min{2NF, No}, (62) 


if the object is either magnetic or dielectric [see (39)]. 


VI. VALIDATION AND CONVERGENCE 


The solution of (25) and (27), which is based on the weak 
form (52) and (53), has been implemented in an extended nu- 
merical code. In this section we consider the validation of the 
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Fig. 5: For testing LEGO-EFIE: surface and volume meshes of two strip- 
dipoles and two dielectric spheres embedded in as many cubic bricks. Data: 
radius of spheres 25 cm, width of dipoles 10 cm, height of dipoles (d) and 
edge of bricks 100 cm, separation of dipoles 220 cm. (After [57].) 


hybrid LEGO-EFIE approach, and we discuss the convergence 
of the solution with the number coupled eigencurrents retained 
for the inversion of (57). 


A. Example of validation 


We examine the simple antenna problem of Fig. 5 in which 
two strip-dipole antennas are symmetrically placed across two 
dielectric spheres (medium ®) from each other. 

The antenna-objects system is immersed in free space 
(medium © and ®© have the same EM parameters), in ac- 
cordance with the general setup sketched in Fig. 1. To model 
the problem with LEGO we embed each sphere inside a cubic 
brick; the center of the sphere coincides with that of the sur- 
rounding brick (see the caption of Fig. 5 for geometrical data). 
The number of basis functions on the dipoles, on the surface 
of each brick and in each sphere are N4 = 15 x 2 = 30, 
2Nr = 1800 and No = 1834, respectively. 

To obtain a reference solution we have solved the problem 
with a previous version of the code [61] for the special case 
where the spheres are comprised of an isotropic dielectric 
material with €3 = 2coI and Fs = uol. For the calculation of 
[S11] in this case we have performed the MoM solution of the 
Poggio-Miller-Chang-Harrington-Wu-Tai SIE [14], [15] on the 
surface of the sphere inside a brick [54] by expanding electric 
and magnetic surface current densities with No = 2 x 294 = 
588 surface divergence-conforming vector basis functions [72] 
of the same type as in (38a), (38b) or (40). 

In assembling the system (53) we have obtained the al- 
gebraic LEGO operator |Luggo] through (59) by employ- 
ing Nc = 30 coupled eigencurrents out of NC max = 
min{2NF, No} = 1800. Likewise, we have constructed the 
system given by [61, Eqs. (9), (14)] by using Nc = 30 coupled 
eigencurrents out of Nomax = min{2Nr,No} = 588. 
Therefore, in both cases the reduced inverse scattering operator 
[Soc]? in (58) has the rank Nc Np = 60. Then, from the 
knowledge of [Ja] and [q] we have computed the 2 x 2 
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Fig. 6: Impedance matrix entries for the problem of Fig. 5: real (a) and imaginary (b) parts versus the electric length of the dipoles; (—) old LEGO-EFIE with 


SIE inside bricks [61], (0) LEGO-EFIE (53) and VIE inside bricks [62]. Data: €, = €2 = £0, E3 = 2eoI, p1 = H2 = H3 
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Fig. 7: Normalized radiated electric field versus the elevation angle (0) for 
the problem of Fig. 5; (—) old LEGO-EFIE with SIE inside bricks [61], (0) 
LEGO-EFIE (53) and VIE inside bricks [62]. Data: see caption of Fig. 6; 
each port excited with Vg = 1 V, d/Ao = 0.467. (After [57].) 


antenna impedance matrix for seventeen frequency samples 
in the range f € [120, 200] MHz, and the normalized radiated 
field at f = 140 MHz. 

The entries Z1ı and Zo, of the antenna impedance matrix 
are plotted in Fig. 6 versus the dipole electric length d/Ag. The 
magnitude of the radiated electric field is plotted in Fig. 7 as 
a function of the elevation angle (0) in the half-planes ¢ = 0° 
and ¢ = 90° at f = 140 MHz. Since the results show perfect 
agreement, we can conclude that the LEGO-EFIE approach 
with VIE is not only viable, but that it has been implemented 
correctly. 


B. Convergence of EEM with No 


It is interesting to investigate how well the solution of (52), 
(53) converges when the number Nc < min{2Npr, No} of 
coupled eigencurrents is increased, because the latter have 
a direct impact on the accuracy of the reduced scattering 
operator [cc] and hence, on [Lygco] in (59). 

Towards this end we have repeatedly solved the antenna 
problem of Fig. 5 at f = 140 MHz (d/o = 0.47) by 
using Nc € {5,10,15, 20, 25,30} coupled eigencurrents. 





uo, No = 30. (After [57].) 
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Fig. 8: Convergence of EEM: L2-norm of EAP k = 1,2, versus number 
of coupled eigencurrents (Vc) for the antenna problem of Fig. 5. Data: see 
caption of Fig. 9. (After [57].) 


Accordingly, we have determined the impedance matrix from 
[Ja] and the radiated field from both [Ja] and [q5]. The L2- 
norms of [q},], k = 1, 2 are shown in Fig. 8, whereas the entries 
411 and Z2; of the impedance matrix are plotted in Fig. 9. 

It is apparent that, for the problem at hand, convergence is 
reached when we employ just No = 15 coupled eigencurrents 
out of 1800. As can be expected, the convergence of the 
radiated fields is even more rapid. This can be ascertained from 
Fig. 10, in which the magnitude of the normalized electric field 
obtained with Nc € {10,20} at f = 140 MHz is compared 
to the same quantity computed with Nc = 30, the latter being 
assumed as a reference solution in view of Fig. 8. 

In Section V we have argued that the coupled (uncoupled) 
eigencurrents of the algebraic scattering operator are asso- 
ciated with the larger (smaller) eigenvalues a9 of [Skk]. 
To support this speculation, in Fig. 11 we have plotted the 
spectrum of [S11] pertinent to the antenna-objects problem of 
Fig. 5 again at f = 140 MHz. As can be seen, the eigenvalues 
decay very fast (exponentially) until they drop down to the 
threshold of numerical noise for double-precision floating- 
point calculations. 

In general, this trend is common to all EM bricks containing 
a PEC [66] or a penetrable object [62], [66]. However, 
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Fig. 9: Convergence of EEM: impedance-matrix entries (Z;1) versus number of coupled eigencurrents (Nc) for the antenna problem of Fig. 5; the markers 





(0) are joined with a dotted line for visualization’s sake. Data: €1 = €2 = €0, €3 = (48% 4 
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Fig. 10: Convergence of EEM: normalized radiated electric field versus 
elevation angle (0) for different numbers of coupled eigencurrents (Nco) for 
the problem of Fig. 5; (0) No = 10, (0) Ne = 20, (—) Ne = 30. Inset: 
close-up of the electric field. Data: see caption of Fig. 9; each port excited 
with Vg = 1 V, d/Xo = 0.467. (After [57].) 














[vs depends substantially on the relative size of brick 
and embedded object and the constitutive parameters thereof. 
Besides, a sharp threshold separating the coupled eigencurrents 
from the uncoupled ones may appear (cf. [66, Figs. 2-4]), if 
2No < 2Npr or No < 2Nr, that is, the number of basis 
functions employed to solve (31) is smaller than the number 
of expansion functions introduced over OD; Actually, in this 
instance the algebraic scattering operator (44) is obviously 
rank-deficient. 


VII. APPLICATION EXAMPLE: PLASMA ANTENNA ARRAY 


A Gaseous Plasma Antenna (GPA) is a plasma discharge 
(usually, in the form of a dielectric tube filled with gas) in 
which the generated plasma behaves as a conducting medium 
when the tube is energized, and this enables the tube to 
transmit and receive EM waves [6], [74]. GPAs have many 
potential advantages over conventional metallic antennas [75], 
namely: 


e They are reconfigurable with respect to radiation pattern, 
frequency, bandwidth [5], [6]. 


29 + 2Z)eo, u1 = u2 = u3 = po. (After [57].) 
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Fig. 11: Convergence of EEM: spectrum of eigenvalues of [S11] as a function 
of their index q = 1,...,2Np for the antenna-objects system of Fig. 5 at 
f = 140 MHz. Inset: close-up of the eigenvalues pertinent to the strongly 
coupled eigencurrents. 


e They can be reconfigured electrically — rather than 
mechanically — on time scales that are on the order of 
microseconds to milliseconds. 

e They have lower thermal noise and minimise signal 
degradation, as they are activated only while communi- 
cation takes place. 

e They are virtually “transparent” above the plasma fre- 
quency, and become “invisible” once turned off. 


Also, GPAs operating at different frequencies do not interfere 
with one another, a property which makes it possible to stack 
arrays of plasma antennas. 

A Plasma Antenna Array (PAA) is a cluster of GPAs and 
possibly conventional metallic elements [74]. In addition to 
having all the advantages of baseline GPAs, PAAs allow (i) 
steering the beam and (ii) improving directivity by adding 
nulls to the radiation pattern. Moreover, an array of plasma 
discharges surrounding a radiating element — which can be 
either metallic or a GPA — effectively realizes a “plasma 
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Fig. 13: Modelling a reconfigurable PAA with LEGO: 3-D triangular tessel- 
lation of two PEC center-fed cylindrical dipoles. Data: height of dipoles (d) 
15 mm, diameter of dipoles 1 mm, separation of dipoles along the x-axis 15 


mm, NA 
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Fig. 14: Modelling a reconfigurable PAA with LEGO: (a) tetrahedral mesh 
of a plasma tube; (b) surface mesh of the brick that embeds the tube. Data: 


diameter of tube 10 mm 


of brick 71 mm, 


20 plasma tubes; 
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blanket” of sorts that can be used to shield unwanted incoming 
EM waves and open “windows” in selected positions so as to 


let outward EM waves through. 
plasma tubes arranged in a circular fashion around the dipoles 


of two PEC center-fed cylindrical dipoles (Fig. 13), whereas 
(Fig. 12a). This PAA is reconfigurable primarily in the sense 


Data: radius of the plasma screen 5 cm. (Also see caption of Figs. 13, 14). 
the plasma blanket is comprised of Np = 


(b) same as before, but with three plasma tubes switched off (Np = 17). 


cylindrical dipoles surrounded by a circular screen of Np 


width of brick 11 mm, height 


height of tube 60 mm, 


> 


= 3056 


, 2Np = 1656. 


No 


that the plasma tubes can selectively be switched on and off 
as needed to create windows for beam-forming and beam- 


steering. For instance, the PAA operating with only Np 


tubes switched on is shown in Fig. 12b. 
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distributed on a circle with radius of 50 mm; each tube is 


(75b) in Appendix B. The cylindrical plasma tubes are evenly 
aligned with the z-axis, 


and (31) takes on the special form given by (74) 


height of 60 mm. The plasma is assumed cold 
collisional, and magnetized. Then, 


coordinates with magneto 


£o. However, since the plasma (medium 


The background and host media are free space; 
E2 


In order to apply the hybrid LEGO-EFIE method, we embed 
each plasma tube in a cuboidal EM brick, as shown in Fig. 14. 
The width and height of the bricks are 11 mm and 71 mm, 
@®) is modelled with a dyadic permittivity, then m = pol, 


respectively. 
therefore £1 
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permittivity €3 reads 











S jD 0 
=e |-jD S Ol, (63) 
0 0 P 
where 
welw — jve) 
s=1-) Bo ee (64a) 
z wilo — jue)” — wee 
yee wre 
D= £ — ; (64b) 
E YH (w = jve)? — w? 
2 
W 
pe (64c) 


al we 
are the Stix parameters [76], with wpe = [neqz/(eome)]*/? the 
plasma frequency, wee = oeqeBo/me the gyrofrequency, ve 
the collision frequency, and og = + the particle charge sign; 
the subscript ‘£’ refers to the plasma species. By letting €3 be a 
function of position, profiles of plasma density, magnetic field, 
electron temperature, and neutral pressure can be handled; 
however, in the numerical experiments described below €3 is 
considered constant. 

The PEC dipoles both have a diameter of 1 mm and a height 
of 15 mm, they are aligned with the z-axis, and are placed at 
a distance of 15 mm along the z-axis (see Fig. 13). They are 
excited in phase with Ve = 1 V [cf. (14)] at the frequency 
f = 10 GHz. 

The surface (volume) basis functions on the boundary of 
a brick (in a plasma tube) are 2Ne = 1656 (No = 3056). 
The current density induced over the dipoles is represented by 
means of N4 = 312 surface basis functions. The problem in 
Fig. 12a is modelled with Np = 20 EM bricks, whereas only 
Np = 17 bricks are employed for the PAA with a window in 
Fig. 12b. 

We have investigated the role played by plasma discharge 
parameters in the two configurations of Fig. 12. Specifically, 
the plasma is made of argon ions (Ar*) and electrons with 
the same uniform density no = 101° m~’. The chosen 
neutral pressure and electron temperature are such that the 
electron collision frequency is ve = 89 GHz, and the collision 
frequency of the argon ions is negligible with respect to ve. 
Two scenarios in terms of magneto-static field have been 
considered: (i) Bo = 0 T (non-magnetized case); and (ii) 
Bo = 0.1 T (magnetized case). With these positions, the Stix 
parameters take on the values listed in Table I. 

The antenna gain function g(@,¢@) relative to the ideal 
isotropic radiator [68] is plotted in Figs. 15 and 16 for the two 
scenarios. For the sake of reference, the gain of the two-dipole 
array in the case when all the plasma tubes are switched off 
(solid blue line) is superimposed to each plot to help appreciate 
the combined effect of the plasma blanket (window) and the 
plasma discharge parameters. 

In particular, from Fig. 15a we notice that the PAA with 
Np = 20 energized plasma tubes exhibits six narrow main 
lobes in the E-plane for œ € {90°, 270°}; correspondingly, 





TABLE I: Plasma parameters for the PAAs of Fig. 12 








Bo = 0 [T] 
S  —7.061727403150083 — j0.013250915072170 
D 0 
P — ~7.061727403150083 — j0.013250915072170 
Bo = 0.1 [T] 
S  —~7.747114953842780 — j0.016822229467298 
D _0.008733555419741 — j2.448490150087137 
P —7.061727403150083 — j0.013250915072170 





the maximum gain is lower that that of the two-dipole array 
in free space. As expected, the plasma window renders the 
gain asymmetric, and enhances the three main lobes in the 
@ = 90° half-plane. In the H-plane (Fig. 15b), the angular 
distribution of power of the PAA with Np = 20 energized 
plasma tubes is similar to that of the two-dipole array in free 
space, but the gain is significantly smaller. When a window is 
opened in the blanket, the beam-forming effect of the active 
plasma tubes becomes evident: A narrow main lobe appears 
in the broadside direction 0 = ¢ = 90°, while the backward 
radiation is reduced. 

Quite similar observations apply to the gain in the second 
scenario (i.e., magnetized plasma tubes) reported in Figs. 16a 
and 16b. In fact, a non-null magneto-static field makes the 
plasma non-reciprocal [cf. (63)], which in turn results in 
asymmetric gain for 6 € {90°,270°} in the E-plane for the 
case Np = 20 plasma tubes. Furthermore, in the H-plane 
the magneto-static field reduces the maximum gain in the 
case of the windowed PAA as compared to the corresponding 
configuration with By = 0 T (Fig. 15b). 

We conclude this discussion with a comment on the effi- 
ciency. While it is clear from Figs. 15 and 16 that the PAA 
under investigation is capable of beam-forming and beam- 
steering, we also note that the maximum gain is lower than 
that of the two-dipole array alone in the H-plane and for some 
directions in the E-plane, owing to the lossy nature of the 
plasma. Therefore, at least in our example higher input power 
should be provided to exploit the reconfigurability afforded by 
PAAs while preserving the level of the maximum gain. 

As regards the computational performance of the hybrid 
LEGO-EFIE approach, the size of the algebraic total inverse 
scattering operator [s}"* in (48) is 2Ne Np = 33120 for 
the PAA with Np = 20 tubes and 2NrNp = 28152 
for the PAA with a window (Np = 17). However, by 
applying the EEM with Nc = 50 coupled eigencurrents out 
of Nomax = min{2Nr,No} = 1656, the calculation of 
[Lueco] in (53) boils down to inverting a matrix [Soc]! with 
a rank NcNp € {1000, 850}. More importantly, the size of 
the system matrix in (53) is just N4 = 312, i.e., determined 
by the number of basis functions used to expand J,4 on 
the dipoles. Furthermore, as all the plasma tubes have the 
same shape and EM properties, only one algebraic scattering 
operator, say [S11], must be computed; therefore, the solution 
of the VIE (31), which is time-consuming, has to be carried 
out only once for each of the scenarios described above. In 
contrast, the direct solution of the problems of Fig. 12 with a 
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Fig. 15: Modelling a reconfigurable PAA with LEGO: antenna gain (natural 
units) for various configurations; (a) E-plane (@ € {90°, 270°}), (b) H-plane 
(6 = 90°). Data: each antenna port excited with Vg = 1 V at f = 10 GHz, 
plasma density no = 1019 m°, plasma magnetizing field Bo = 0 T. 


VSIE formulation (e.g., [7]) would require the solution of a 
system with rank NoNp + Na € {61432, 52264}. 


VIII. CONCLUSION 


The hybrid LEGO-EFIE method uses EM bricks in tan- 
dem with the standard EFIE to model complex problems, 
which involve conventional PEC antennas or scatterers and 
media with anisotropic constitutive parameters. The strategy 
is particularly well-suited when the problem at hand comes 
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Fig. 16: Modelling a reconfigurable PAA with LEGO: antenna gain (natural 
units) for various configurations; (a) E-plane (ġ € {90° , 270° }), (b) H-plane 
(8 = 90°). Data: each antenna port excited with Vg = 1 V at f = 10 GHz, 
plasma density no = 1019 m?, plasma magnetizing field Bo = 0.1 T. 


naturally decomposed in many separated parts. The calculation 
of the scattering operators of the bricks with the baseline MoM 
entails the inversion of one small VIE at a time, which is 
extremely competitive over the solution of coupled VSIEs for 
the entire problem. Besides, the very introduction of scattering 
operators enables us to reduce the degrees of freedom, because 
expanding the unknown equivalent current densities on the 
boundary of a brick usually requires fewer basis functions than 
are necessary to solve the VIE for the embedded anisotropic 
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object. As the modularity of the approach enables us to 
combine and recycle the EM bricks with few restrictions, the 
time required to analyze similar structures (such as the PAAs 
addressed in Section VID) can also be considerably reduced. Fi- 
nally, the usage of the eigencurrents of the scattering operators 
as locally entire-domain basis functions over the boundaries of 
the bricks enables us to handle and compress the weak form 
of the LEGO functional equations to derive the modified EFIE 
efficiently. 


APPENDIX A 
PROPAGATORS USED IN LEGO 


Listed in this appendix are the propagators involved in the 
definition of Skk, Tkn, (Toa), and (Pao). 

In what follows, kı denotes the wavenumber in medium ©, 
r (r’) denotes the field (source) point, I, is the transverse 
unit dyadic tangential to OD;,, or Sa. The surface ‘del’ 
operator is defined as V, = I, - V, and Vi = —V, means 
differentiation with respect to r’. The specification ‘PV’ before 
an integral sign means that integration is performed in the 
Cauchy principle value sense, while d?r’, dèr’ represent the 
area and the volume elements, respectively. The unit normal 
hy, to OD; points inward Dy (Fig. 1b). 


A. Brick-to-brick interaction 
The propagators Pin Pèp Pkn are 2 x 2 abstract matrices 
whose entries are dyadic integro-differential operators involv- 
ing the Green’s function G(R) in (17). 
Pir H Dk => Dk 
(from incident currents to incident fields) 


(Pip) {0} = -jki |  PrGR)T {0} 








a d?r’Gi(R)Vi- {0}, re ODz, (65a) 
OD; 
(Pix)i2{o} = —PV pat) x Is- {o} 
k 
1 A 
+ zô xI,-{o}, reOD,, (65b) 
Pidato} =v far’ V.G(R) x Ts- {0} 
k 
1 = 
— zô x Is-{o}, reððDp, (65c) 
(Pix)a2{o} = —(Pix)ar {>}, (65d) 
Pry : OD; —> OD; 
(from scattered currents to scattered fields) 
(Pé,)ar{o} h I rG RT,- {0} 
k 
IVs f a2r'GI(R)V}- {0}, r€ODg, (66a) 
> ky ap r Gi s OF, r k; a 
k 


(Pg, )i2{o} = -Pv d?r’ViGi(R) x I, - {0} 


OD 
1 = 
— aie xI,-{o}, reOD,, (66b) 
(Pèp)21{0} = evf d?r'V,G1(R) x I, . {o} 
OD, 
1 = 
+ ik xI,-{o}, reODz, (66c) 
(Pr)o2to} = — (Pře) {0}, (66d) 
Prk x OD, =} OD, n # k 
(from scattered currents to incident fields) 
(Pin )11 {0} = (Pkn)i12{0} = O- {0}, (67a) 
(Pinlarto} sev d?r'V.Gx(R) x T- {0} 
ODn 
1 = 
+ zô xI,-{o}, reOD,, (67b) 


(Pin )a2{o} = iaf d?r’Gi(R)I, - {0} 





+f d?r’Gi(R)Vi-{o}, re ODz, (67c) 
ki Jap, 


The entries (Pkn)i1 and (Pkn)ı2 are null, because the 
propagator Pi, [cf. (36a)] requires only the magnetic field 
to determine the equivalent incident currents. 


B. Brick-to-object and object-to-brick interaction 


In the special case when the anisotropic body is either 
dielectric (with g, = pD or magnetic (with €; = J), 
the propagators Pox (Pko) reduce to 1 x 2 (2 x 1) abstract 
matrices whose entries are dyadic integro-differential operators 
involving the Green’s function G(R) in (17). The contrast 
dyadics @,.(r) and @),(r) are defined in (30a) and (30b). 


Anisotropic dielectric body 
Pok $ OD; =7 Vo 
(from incident currents to incident electric flux density) 


(Pox) fo} = -jki i) Ëra) {0} 





= w. ; d?°r'Gi(R)V;- {0}, re Vo, (68a) 
Dk 
(Por)ia{o} = -Pv d?r'VGi(R) xT, {0} 
OD, 
1 = 
= gtr xI,-fo}, reVo, (68b) 


Anisotropic magnetic body 
Por : OD, =} Vo 
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(from incident currents to incident magnetic flux density) 
(Pox)21 = -vf d?r’ViGi(R) x Is- {0} 
OD, 


1 = 
-zô xT. {0}, re Vo, (69a) 


(Pox )22{0} = jki d?r’G(R)Is - {0} 
Dk 


Vs f a?r Gi(R)V!- fo}, 
kı Jap, 





+ r€ Vo, (69b) 
Anisotropic dielectric body 
Pko í Vo => Dk 


(from total electric flux density to scattered fields) 


(Pko)u{o} = 0: {0}, (70a) 


(Pro)z1{0} = jki f d?r’ Te - VGi(R) x Telr’) - {0}, 


Vo 


reðD,, (10b) 


Anisotropic magnetic body 
Pro: Vo > OD; 
(from total magnetic flux density to scattered fields) 


(Pro)i2{o} = 0- {o}, (71a) 
(Pko)22{0} = k2 i d?r’Gi(R) Is -œa (r) - fo} 
+V, | d3r’Gi(R)V' - [@n(r’)-{o}], re OD, (11b) 


Vo 


The entries (Px.o)11 and (Pk o)12 are null, since the propaga- 
tor Pè, [cf. (34)] requires only the magnetic field to determine 
the equivalent scattered currents. 


C. Antenna-to-brick and brick-to-antenna interaction 


The propagators Pa (Pax) are 2 x 1 (1 x 2) abstract ma- 
trices, whose entries are dyadic integro-differential operators 
involving the Green’s function G(R) in (17). 


Pra SA =% OD, 
(from antenna current to incident fields) 


(Pra )is{o} = 0- {0}, 


(Pra )ai{o} = p d?r’V.Gi(R)xI,-{o}, 


(72a) 
r € OD, (72b) 


Par: OD, > SA 
(from scattered currents to incident fields) 


d?r’G,(R) I, - {o} 
OD, 


(Pax)ii{o} = —jki 





~Ms f 2a (RV! - fo}, reSa, (73a) 
ki Javx 
(Par)i2{o} = d?r'V,Gi(R)xI,-{o}, re Sa, (73b) 


OD; 


The entry (Pka)ıı is null, because the propagator Pi, 
in (36a) requires only the magnetic field to determine the 
equivalent incident currents. Provided the antenna and the 
conducting parts do not touch the bricks, it is not necessary 
to take the principal value of the integrals (72b) and (73b) nor 
do they contribute a residue, for the kernel is never singular. 


APPENDIX B 
EXPLICIT FORMS OF (31) 


In the special instance of either dielectric (with #3 = 11) 
or magnetic (with €3 = aD anisotropic object, Xoo in (31) 
reduces to a volume integro-differential operator involving the 
Green’s function G;(R) in (17) and either @(r) or @),(r). 
Also, go and Fi become one-element abstract column vectors. 


Anisotropic dielectric body 
Xoo : Vo > Vo 
(from flux density to flux density) 





Xoo{o} = 22] E k? $ d?r’Gı (R) Telr’) - {0} 
-V ; d3r'Gi(R)V' - [a&e(r’)- {0}, reVo, (74) 
do = Do/(Er. Vm); (75a) 
Fi = Di/(e1vm); (15b) 


Anisotropic magnetic body 
Xoo : Vo > Vo 
(from flux density to flux density) 





Xeo(0} = ES nese [ Eroa) o 


-V | d’r'Gi(R)V'-[@a(r’)-{o}], reVo, (76) 
Vo 

qo = Bom / H1, a) 

FÌ = Bİ Vin/m, (7b) 
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